今回紹介する手法は、油で飽和された岩石を水で掃攻する際に飽和率の時間変化を求める事が出来るというものである。
私はあまり理論の歴史に詳しくはないが、書籍などによるともとは1942年にBuckleyとLeverettという人たちが提唱した手法とのこと。 もっと複雑な状況も扱えるようだが、勉強不足のため今回は最も単純な状況の理論をまとめる。 また、個人的に納得しきれていないところもあるため、今後改良しようと考えている。
水平方向（$x$方向）のみに流体が流動し、毛細管圧力の影響はないと仮定する。ダルシーの法則より水の流束$u_w$と油の流束$u_o$は
$$u_w=-\frac{kk_{rw}}{\mu_w} \frac{\partial p}{\partial x}$$$$u_o=-\frac{kk_{ro}}{\mu_o} \frac{\partial p}{\partial x}$$と表せる。ここで、$k$は絶対浸透率、$k_{rw}$は水の相対浸透率、$k_{ro}$は油の相対浸透率、$\mu_w$は水の粘性係数、$\mu_o$は油の粘性係数、$p$は圧力である。
水のフラクショナルフロー$f_w$は単純に全体の流束のなかで水の流束が占める割合であり、以下のように定義される。
$$f_w = \frac{u_w}{u_w+u_o}$$これは、ダルシーの法則を用いると
$$f_w=\frac{\frac{k_{rw}}{\mu_w}}{\frac{k_{rw}}{\mu_w}+\frac{k_{ro}}{\mu_o}} \tag{1}$$と書き換えられる。
ここで質量保存の法則より、
$$\begin{align} \frac{\partial}{\partial x} (\rho_w u_w)+\frac{\partial}{\partial t}(\rho_w \phi S_w) &=0 \\ \frac{\partial}{\partial x} (\rho_o u_o)+\frac{\partial}{\partial t}(\rho_o \phi S_o)&=0 \end{align}$$$\rho_w$、$\rho_o$、$\phi$は一定と仮定すると
$$\begin{align} \frac{\partial u_w}{\partial x} +\phi\frac{\partial S_w}{\partial t}&=0 \tag{2}\\ \frac{\partial u_o}{\partial x} +\phi\frac{\partial S_o}{\partial t}&=0 \tag{3} \end{align}$$辺々加えると
$$\frac{\partial }{\partial x}(u_w + u_o) +\phi\frac{\partial}{\partial t}(S_w + S_o)=0$$$S_w + S_o=1$なので
$$\frac{\partial}{\partial x} (u_w + u_o)=0$$が成立する。非圧縮性なので当たり前であるが、$u_w$と$u_o$の和は位置によらず一定となる。$u_w$と$u_o$の和を$u_T$とし、$u_T$は時間に関して一定であると仮定する。 式（2）の水に関する質量保存の法則は$f_w$を用いて以下のように書き換えられる。
$$u_T \frac{\partial f_w}{\partial x}+\phi \frac{\partial S_w}{\partial t} =0$$式（1）のように$f_w$は相対浸透率と粘性係数で表せた。粘性係数は一定であるとすると相対浸透率は飽和率のみの関数として表せる。 飽和率は$S_w+S_o=1$という条件を満たすため、$f_w$は$S_w$の一変数関数とみなす事が出来る。このことから、
$$\frac{u_T}{\phi} \frac{df_w}{dS_w}\frac{\partial S_w}{\partial x}+ \frac{\partial S_w}{\partial t} =0 \tag{4}$$と変形できる。
今$S_w$は$x$と$t$の関数であるが、解きやすいように変数変換を行う（特性曲線法）。 $S_w$を$s$と$\tau$の関数とし、$S_w$を$s$で微分する。
$$\frac{\partial S_w}{\partial s}=\frac{\partial S_w}{\partial x}\frac{\partial x}{\partial s}+\frac{\partial S_w}{\partial t}\frac{\partial t}{\partial s} \tag{5}$$次に、式（5）が式（4）と同じ形になるように以下のように$s$と$\tau$を定めていく。これ以降$g_i(\tau)(i \in \mathbb{N})$は何らかの$\tau$の関数を表す。
$$\begin{align} \frac{\partial S_w}{\partial s}&=0 \tag{6}\\ \frac{\partial x}{\partial s}&=\frac{u_T}{\phi} \frac{df_w}{dS_w}\tag{7}\\ \frac{\partial t}{\partial s}&=1 \tag{8} \end{align}$$最初に式（6）は$S_w$は$s$の関数ではないことを表しているので、
$$S_w=g_1(\tau)$$式（7）の右辺は$S_w$のみの関数であるが、$S_w$は$s$の関数ではないので、右辺自体が$s$の関数ではない。よって
$$x=\frac{u_T}{\phi} \frac{df_w}{dS_w}s+g_2(\tau)$$が得られる。式（8）は
$$t=s+g_3(\tau)$$となるはずである。$s=0$の時$t=0$、$x=\tau$となるようにすると
$$\begin{align} g_2(\tau)&=\tau\\ g_3(\tau)&=0 \end{align}$$となる。以上をまとめると
$$\begin{align} S_w&=g_1(\tau) \tag{9}\\ x&=\frac{u_T}{\phi} \frac{df_w}{dS_w}s+\tau \tag{10}\\ t&=s \tag{11} \end{align}$$となり、式（10）と式（11）より
$$\tau=x-\frac{u_T}{\phi} \frac{df_w}{dS_w}t \tag{12}$$である。後程触れるが共通の$x$、$t$に対して異なる$S_w$で同じ$\tau$を示す可能性がある。実は最終的にこの解が適用されるのは$S_w$と$\tau$が一対一に対応する場合に限定されるため、ここでは$S_w$と$\tau$が一対一に対応すると仮定する。（正直このあたりは私は完全に理解できていない。後にショックフロントを導入するのだが、「微分方程式に微分できない解が存在すること」に関してはまだ勉強中であるため、ごまかしている。） 式（9）と式（12）が同時に満たされるように$g_1$ を選ぶ必要があるため、$x$、$t$が定数であれば式（9）と式（12）は等価な式である。式（12）はある$S_w$、$x_1$、$t_1$に対して
$$\tau=x_1-\frac{u_T}{\phi} \frac{df_w}{dS_w}t_1 \tag{13}$$とした時、任意の$t_2( \neq t_1)$に対しも式（13）が成立するような$x_2$が存在するはずである。
$$\tau=x_2-\frac{u_T}{\phi} \frac{df_w}{dS_w}t_2 \tag{14}$$式（13）と式（14）を辺々引き、整理すると
$$\frac{x_2-x_1}{t_2-t_1}=\frac{u_T}{\phi} \frac{df_w}{dS_w} \tag{15}$$となる。式（15）の右辺は$S_w$のみの関数であり、$S_w$が一定となるとき$x_2-x_1$と$t_2-t_1$の比は一定になることを示している。 つまり、ある$S_w$は一定の速度で進行していくことを示している。
これまでの理論では式（12）の$x$と$t$を固定した時に$S_w$と$\tau$が一対一の対応をしていなければならなかった。ここで、一般的な$k_{rw}$、$k_{ro}$、$f_w$、$\displaystyle \frac{df_w}{dS_w}$のグラフを以下に示す。
式（15）より、ある$S_w$の進行速度は$\displaystyle \frac{df_w}{dS_w}$に比例することが分かるが、$\displaystyle \frac{df_w}{dS_w}$は異なる$S_w$でも同じ値を示す可能性がある。つまり、異なる$S_w$が同じ速度で進む可能性が出てきてしまう。これは物理的におかしな現象である。 前にも述べたように限られた範囲でしかこの$\displaystyle \frac{df_w}{dS_w}$と$S_w$の関係は使用しない。 これ以降どこまでの範囲を使用するかに関して議論を進める。
初期の水飽和率を$S_{wi}$とし、岩石を時刻$t=0$、位置$x=0$から流束$u_T$で水で掃攻することを想定する。また、岩石の断面積を$A$、残留油飽和率を$S_{or}$とする。 $x=0$から水を圧入しているので$x=0$では瞬時に水飽和率が$1-S_{or}$になるはずである。 圧入した水が到達している先端の位置（ショックフロントと呼ばれる）を$x_f$、水飽和率を$S_{wf}$とし、$x=0$から$x=x_f$までの平均の水飽和率を$\overline{S_w}$とする。非圧縮性流体を仮定しているため、圧入した水の分だけ水飽和率が増えるので
$$\begin{align} u_TtA &=A \phi \int_0^{x_f}(S_w-S_{wi}) \ dx \\ &=A \phi \int_0^{x_f} S_w \ dx - A \phi S_{wi} x_f \\ &=A \phi \overline{S_w} x_f - A \phi S_{wi} x_f \\ &=A \phi x_f (\overline{S_w} - S_{wi}) \end{align}$$もう少し変形をすると
$$ \frac{\phi}{u_T}\frac{x_f}{t}=\frac{1}{\overline{S_w} - S_{wi}} \tag{16} $$となり、$\displaystyle \frac{x_f}{t}$は水が到達している先端の速度なので式（15）より
$$ \left( \frac{df_w}{dS_w} \right)_{S_w=S_{wf}}=\frac{1}{\overline{S_w} - S_{wi}} \tag{17} $$となる。次に先端位置における時間$\Delta t$での物質収支を考える。
$$ Au_T(f_w(S_{wf})-f_w(S_{wi}))\Delta t = A\phi (S_{wf}-S_{wi}) v_{sf}\Delta t $$ここで$v_{sf}$は先端（ショックフロント）の速度である。式（15）を用いると、
$$ \left( \frac{df_w}{dS_w} \right)_{S_w=S_{wf}}= \frac{f_w(S_{wf})-f_w(S_{wi})}{S_{wf} - S_{wi}} \tag{18} $$が得られる。
式（17）と式（18）は非常に興味深い。これらを用いる事でショックフロントの水飽和率$S_{wf}$とショックフロントまでの平均の水飽和率$\overline{S_w}$を求める事ができる。 具体的に言うと横軸$S_{w}$縦軸$f_w$のグラフで$(S_{wi}, f_w(S_{wi}))$から接線を引く。接点がショックフロントの水飽和率となる。以下のグラフの例では、初期の水飽和率が不動水飽和率の時を想定している。
これまでの理論よりショックフロントの水飽和率を算出する事が出来た。あとは、ショックフロントまでの水飽和率の分布を$\displaystyle \frac{df_w}{dS_w}$に従って求めればよい。
上では述べなかったが、実はショックフロントの水飽和率は以下のAとB領域の面積が等しくなるようになっている。
これは、単純に以下のように示せる。
$$\begin{align} \int_{S_{wi}}^{S_{wf}}\left( \frac{df_w}{dS_w} - \left( \frac{df_w}{dS_w} \right)_{S_w=S_{wf}} \right) \ dS_w &= \int_{S_{wi}}^{S_{wf}}\frac{df_w}{dS_w} \ dS_w - \left( \frac{df_w}{dS_w} \right)_{S_w=S_{wf}} (S_{wf}-S_{wi}) \\ &= f_w(S_{wf})-f_w(S_{wi})- \left( \frac{df_w}{dS_w} \right)_{S_w=S_{wf}} (S_{wf}-S_{wi}) \\ &= 0 \end{align}$$最後は式（18）の関係式を利用している。
また、ショックフロントの位置までの任意の位置における平均水飽和率もグラフを用いて求める事ができる。 少々技巧的なため後回しにしたが、初めに紹介した物質収支を考えた方法はこれから示す方法の特殊な状況である。
水飽和率が$S_{w1}$の位置を$x_1$とし、$x_1$までの平均の水飽和率を$\overline{S_{w1}}$、$x_1$における水のフラクショナルフローを$f_{w1}$とする。
$$\begin{align} \overline{S_{w1}} &= \frac{1}{A \phi x_1}\int_0^{x_1}A \phi S_w \ dx \\ &= \frac{1}{x_1}\int_0^{x_1} S_w \ dx \\ &= \frac{1}{x_1}\left( \left[ xS_w \right]_0^{x_1}- \int_0^{x_1} x \frac{\partial S_w}{\partial x} \ dx \right) \quad \mbox{（部分積分を使用した）}\\ &= \frac{1}{x_1}\left( x_1S_{w1} - \int_0^{x_1} x \frac{dS_w}{df_w}\frac{\partial f_w}{\partial x} \ dx \right) \\ &= S_{w1} - \frac{1}{x_1}\int_0^{x_1} x \frac{1}{\frac{df_w}{dS_w}}\frac{\partial f_w}{\partial x} \ dx \\ &= S_{w1} - \frac{1}{x_1} \int_0^{x_1} x \frac{u_T t}{\phi x}\frac{\partial f_w}{\partial x} \ dx \quad \mbox{（式（15）を使用した）} \\ &= S_{w1} - \frac{u_T t}{\phi x_1} \int_0^{x_1} \frac{\partial f_w}{\partial x} \ dx \\ &= S_{w1} - \frac{1}{\left(\frac{df_w}{dS_w}\right) _{S_w=S_{w1}}} \int_{1}^{f_{w1}} df_w \quad \mbox{（式（15）を使用した）}\\ &= S_{w1} - \frac{f_{w1}-1}{\left(\frac{df_w}{dS_w}\right) _{S_w=S_{w1}}} \end{align}$$以上より
$$ \left(\frac{df_w}{dS_w}\right) _{S_w=S_{w1}}=\frac{1-f_{w1}}{\overline{S_{w1}}-S_{w1}} $$